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Abstract 

A family of continuous piecewise linear finite elements for thin plate problems is 
presented. We use standard linear interpolation of the deflection field to reconstruct a 
discontinuous piecewise quadratic deflection field. This allows us to use discontinuous 
Galerkin methods for the Kirchhoff-Love plate equation. Three example reconstructions 
of quadratic functions from linear interpolation triangles are presented: a reconstruction 
using Morley basis functions, a fully quadratic reconstruction, and a more general least 
squares approach to a fully quadratic reconstruction. The Morley reconstruction is shown 
to be equivalent to the Basic Plate Triangle. Given a condition on the reconstruction op¬ 
erator, a priori error estimates are proved in energy norm and U norm. Numerical results 
indicate that the Morley reconstruction/Basic Plate Triangle does not converge on unstruc¬ 
tured meshes while the fully quadratic reconstruction show optimal convergence. 
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1 Introduction 

The Kirchhoff-Eove plate equation is a fourth order partial differential equation modeling the 
deflection of thin plates. To approximate solutions to this equation using standard finite ele¬ 
ment methods finite element spaces are required. The difficulty of creating sueh spaees on 
unstruetured triangulations is a well known problem. A possible element is the conforming 
Argyris triangle Q whieh use a fifth order polynomial approximation. Noneonforming options 
inelude the Morley triangle ffTOl and more reeently diseontinuous Galerkin (dG) methods |]3[H1- 
While it is elear that higher order elements are in many ways superior for modeling the plate 
equation, an advantage of low order elements lies in modeling eomplex domains using few 
degrees of freedom. With the extension to shells and the desired eonformity when eombining 
shells and volumes the advantages of low order elements that only feature displaeement de¬ 
grees of freedom become obvious. While this is a possibility when using dG methods, eurrent 
formulations [[3 [H require at least pieeewise quadratic polynomials to yield accurate results. 
The focus of this paper is accurate modeling of the plate equation using a eontinuous piecewise 
linear defleetion field. 

Several authors have tried to develop finite element methods for thin plate modeling using 
a eontinuous pieeewise linear defleetion field. Sinee most terms in the variational formulation 
then vanish there is a need to diseretely approximate higher order quantities to retain suffleient 
information. Therefore a eommon trait for this elass of elements is that patehes of elements 
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are used for these approximations. Nay and Utku [ITT]| used a patch of elements to reconstruct 
a quadratic deflection field on each element using least squares approximation. Barnes [|2l 
introduced a facet triangular plate element where the normal curvature to each edge is approx¬ 
imated from the change in normal gradient to neighboring elements. In a similar approach 
Hampshire [|3 derived a plate element where the stiffness was represented by torsional springs 
at each edge. Also based on the idea of torsional springs at element edges Phaal and Calladine 
llT4l [TSil presented a family of facet plate and shell elements which use quadratic polynomial 
reconstruction to calibrate the spring coefficients. By using a mixed interpolation technique in 
combination with finite volume concepts Onate and Cervera IfT^ and Onate and Zarate [fT3l 
proposed a procedure for deriving linear thin plate and shell elements. 

In this paper we present a framework for constructing continuous piecewise linear finite 
elements for the Kirchhoff-Love plate equation. The fundamental idea is to use patches of a 
continuous piecewise linear function to reconstruct a discontinuous piecewise quadratic func¬ 
tion which is used in a dG formulation. We apply the framework for reconstructions in a finite 
element formalism presented in [[3]| to a general dG method for the Kirchhoff-Love plate equa¬ 
tion [[HI. Three example reconstructions are presented and related to existing elements. Given 
a condition on the reconstruction operator we prove a priori error estimates in the energy norm 
and in the norm. 

The remainder of this paper is organized as follows; in Section 2 we present the Kirchhoff- 
Love plate model and the discontinuous Galerkin method using piecewise quadratics contin¬ 
uous at the nodes, in Section 3 we present three reconstructions from continuous piecewise 
linears into piecewise quadratics, in Section 4 we prove a priori error estimates, and in Section 
5 we present convergence studies and numerical examples. 


2 The Plate Model and dG Method 


2.1 The Kirchhoff-Love Plate Model 

The Kirchhoff-Love equilibrium equation governing the deflection of a thin elastic plate occu¬ 
pying a plane domain Q takes the form: Given /, find the deflection u such that 


~ f ^ 


( 2 . 1 ) 


where we use the summation convention and the comma sign indicates differentiation. The 
relationship between moments Oij and curvatures Kij is given by 

Gjj = XAu8ij^ i^j= 1^2, ( 2 . 2 ) 


where dij is the Kronecker delta, A is the Laplacian, A and /i are Lame parameters, and Kij are 
curvatures defined by Kij{u) = u^ij. Using Poisson’s ratio v and bending stiffness D we can 
write the Lame parameters A = Dv and /i = D(1 — v). The bending stiffness of the plate is 
defined by 


Ep^ 

12(1-v2) 


where E is Young’s modulus and p is the thickness of the plate. 


(2.3) 
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Let n = (ni,n 2 ) be an outwards unit normal to the boundary F = dQ. and let t = = 

(^ 2 , —n\) be a tangent to F. To define the boundary conditions we need the following quantities 


U^n = UjHj (2.4) 

u,t = Ujtj (2.5) 

^nn — ( 2 - 6 ) 

^nt — ( 2 - 7 ) 

T = OijjHi-]-M„t,t (2.8) 


where and are normal and tangential gradients, and M„t are bending and twisting 
moments, and T is the transversal force. 

We split the boundary into three disjoint parts F = FcUF 5 UFf’ and let these parts define 
a clamped boundary, a simply supported boundary, and a free boundary. Let the set of angular 
comers on F/? be denoted Xp. The boundary conditions read 


u = Uji = 0 

onFc 

(2.9) 

U = Myifi = 0 

onF^ 

( 2 . 10 ) 

Mnn = T=0 

on Ff- 

( 2 . 11 ) 


at Xp 

( 2 . 12 ) 


where and denote the normal and tangent of F at respective sides of an 

angular corner. 

Let H^{co) denote the Sobolev space of order 5 on the set COCQ., with norm || ■ || ^ and semi¬ 
norm (a defined for m < s. Introducing the following function space where the essential 
boundary conditions are imposed 

W = {v G : V = = 0 on Fc, v = 0 on F 5 } (2.13) 

we recall that the standard variational statement reads: Find w G W such that 


{Oij{u), Kij (v)) = (/, v) for all V G W 


(2.14) 


The calculations leading to this variational statement will be performed in Section 2.L2[ albeit 
on an element level. 


2.1.1 The Mesh and Discontinuous Space 

Let 1C = {K} he a triangulation of Q. into geometrically conforming shape regular triangles. We 
denote the diameter of element K by hp and the global mesh size parameter by h = maxp^p hp. 
Further, let the mesh be quasi-uniform such that 

ch<hp<Ch for all ii: (2.15) 

where c and C are mesh independent constants. The set of edges in the mesh is denoted by 
S = {E} and the set of nodes in the mesh is denoted by V = V(/C) = V{S). We split 8 into 
disjoint subsets 


8 = SiVJ 8c U £^5 U 8p 


(2.16) 
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where Si is the set of edges in the interior of Q., Sc is the set of edges on Pc, ete. Further, with 
eaeh edge we assoeiate a fixed unit normal ng and a eorresponding unit tangent tc sueh that for 
edges on the boundary he is the exterior unit normal. On eaeh node dE belonging to edge E 
we define dnc = 1 if ?£■ points outwards from E and dnc = — 1 is points inwards to E. 

For reasons that become evident when we define the reconstruction operators we make a 
special construction: for every exterior edge E G S\Si we add a ghost element outside the 
domain by placing an additional a degree of freedom, a ghost node, such that the ghost element 
becomes anti-symmetric to the interior element, see Figure |l(b)[ We denote the set of ghost 
elements by JCc- 

Next we define a number of function spaces: Let CVi(S) denote the space of continuous 
piecewise linear functions with support on a set of elements S 

CVi{S) = {veC^{a): v|/c G Pi(i^) for all i^G 5} (2.17) 

and let CVi denote the space of continuous piecewise linear functions with support on /C U ICq 
and zero on the clamped and the simply supported boundary 

CPi = {vgCPi(/CU/Cg) : v = 0onxG£’cU^5} (2.18) 

Furthermore, let W 2 denote the space of discontinuous piecewise quadratic polynomials 

VV 2 = {v : v\k G P 2 (P) for all P G /C} (2.19) 


and finally let VW denote the space of discontinuous piecewise quadratic polynomials that 
are continuous at the nodes and zero on nodes associated with the clamped and the simply 
supported boundaries 


VW = {v G VV 2 : V continuous in x G V, v = 0 in x G V{Sc U Ss) } 
To formulate our method we will use the following notation for the average 

(v+ -I-v“)/2 E & Si 


(v) = 


and for the jump 


V = 


EeS\Si 

— v~ E & Si 
^ EeS\Si 


of a function v at an edge E, where = lim£^o+ eng) with x G P. 


( 2 . 20 ) 


( 2 . 21 ) 


( 2 . 22 ) 


2.1.2 Variational Formulation on an Element 


As a motivation for the dG method we will here derive a variational formulation on each ele¬ 
ment. We multiply ( |2.1[ ) by a test function v G El^ = El^{EL) and integrate over K. Applying 
Green’s formula two times gives 


= ~ (2.23) 

= {^iji ^,ij)K ~ {MnmV^n)dK ~ iMnt^ V,t)dK “h dK 


5 







The final publication is available at: link.springer.com 
Numerische Mathematik (2012) 121:6597 
DOI 10.1007/s00211-011-0429-5 


where we use that vj = j + v^ttj in the last equality. 
Partial integration along an edge segment E gives 


(2.24) 


Combining ( |2.23| ) and ( |2.24[ ) we have the following variational formulation on the element 
level 


{<Jij{u),Kij{v))K- ^ \ {.Mnn:V^n)E-iT,v)E + {M„t,VnsE)dEj={f^v)K (2.25) 

EcdK 


for all V G H^. 


2.1.3 Discrete Moments and Corner Forces 


By giving definitions of the bending and twisting moments and the transversal foree on ele¬ 
ment edges for funetions in WV whieh is eonsistent for funetions in we ean extend the 
elementwise variational statement ( |2.25[ ) to a variational sta teme nt on -\- WV. Following 
the procedure in [[8]| and motivated by the proof of Lemma 4.5 below we for v E + WV 


introduce the following definitions of these quantities on each element edge E G S unless pre¬ 
viously defined by boundary conditions: 


^nniy) i^nniy')) 
r(v) = (r(v)) 
Mntiv) = {Mntiv)) 


■ph ^Po[v,„] 


(2.26) 

(2.27) 

(2.28) 


where /3 is a positive parameter and Pq is the projection onto the space of constants. Using 


these definitions in ( |2.25| ) and summing over all elements K E )C yields a variational statement 
on H^ + WV. 

Due to the nodal continuity of El^ + WV terms containing the twisting moment will vanish 
on all interior edges. On the boundary pointwise twisting moments will appear where the 
boundary is not smooth, but given the homogeneous boundary conditions these terms will be 


zero on Fc U r 5 as V = 0, and also zero on F/? due to ( |2.12[ ). 

The resulting variational statement is nonsymmetric but we may symmetrize the variational 
statement without affecting consistency as the added terms become zero for the exact solution. 
Next we present the resulting variational statement on -|- WV. 


2.1.4 Extended Variational Statement 

The extended variational statement reads: Find u E + WV such that 

a{u,v)=l{v) for all vEH'^ + WV (2.29) 
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where the bilinear form is defined by 

aiv,w) = '^{Oij{v),Kij{w))i 


)K 


KeK. 


y{{^nn{v)) ,[w^n])£ + {[v,n] , {Mnn{w)))^ 

EeSMSsUSF) 

- j5{h^^Po[v^n],Po[w,n])E 


+ E (((7’M),M)^ + ([v],(r(w))), 

Ee£\£F 


where j8 is a real parameter and the linear funetional is defined by 

Kv) = (/,v) 

We now move on to formulate the dG method. 


(2.30) 


(2.31) 


2.2 The dG Method with Piecewise Quadratics Continuous at Nodes 

The dG method for the plate equation with pieeewise quadratie funetions eontinuous at the 
nodes can now be formulated as follows: Find U E WV such that 


a{U,v) = l{v) for all V G'D'PV 


(2.32) 


where the bilinear form is given by ( |2.30| ) and the linear functional is given by ( |2.31[ ). Note 
that the last sum in the bilinear form ( |2.30[ ) gives no contribution as {T (v)) =0 for v E WV. 
The boundary condition w „ = 0 on Fc is weakly enforced via the /3 penalty term while the 
condition M = 0 onrcUr 5 is strongly enforced at the nodes. 

For a more general dG method for the plate equation without the restriction to nodal conti¬ 
nuity and piecewise quadratics in the approximation of the deflection field we refer to |l8l . 


2.3 The dG Method with Embedded Continuous Piecewise Linears 

To formulate our method using a continuous piecewise linear deflection field we use the frame¬ 
work presented in [|3l for using reconstructions in a finite element formalism. We let TZ be 
a reconstruction operator which embeds the space of continuous piecewise linear polynomial 
functions CVi into the space WV of discontinuous piecewise quadratic polynomials continu¬ 
ous at the nodes: 


n-.CVi^WV (2.33) 

Also let the following criterion on the reconstruction operator hold: For v E CV\ 

V = 77.V, for all x G V (2.34) 

The discontinuous Galerkin method with embedded continuous piecewise linear functions 
takes the following form: Find U G CV\ such that 

a{nu,nv) = /(7^v), for all v G CVi (2.35) 
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where a(-,-) and /(■) are defined in ( |2.30[ ) and ( |2.31[ ). The clamped boundary condition is 
weakly enforced by the /3 penalty parameter on TZU. As U coincides with TZU at the nodes we 
can choose to strongly enforce M = 0 onrcUr 5 directly on U. 


3 Examples of Reconstruction Operators 

In this section we consider three reconstruction operators in the presented framework, all of 
which embed continuous piecewise linear functions into WV. To reconstruct a quadratic 
function on an element K these operators use the vertex information in a patch of elements. 
In the first example we reconstruct into the space of the quadratic Morley basis functions, 
which is the subspace of functions in WV that have continuous normal derivative at element 
edge midpoints. We show that this method is equivalent to the Basic Plate Triangle presented in 
iiEiini. The second example is a fully quadratic reconstruction into WV using a four element 
patch. In the last example reconstruction we handle special cases where the fully quadratic 
reconstruction breaks down due to the mesh configuration. A least squares approach to fully 
quadratic reconstruction is used to allow larger patches when fully quadratic reconstruction 
from a four element patch fails. 

3.1 Patch of Elements 

To reconstruct a complete quadratic polynomial six independent degrees of freedom are re¬ 
quired. Thus, a patch of continuous piecewise linear elements is needed to represent sufficient 
information. We denote the patch that is used for reconstructing a quadratic function on ele¬ 
ment K by JV{K) and let it consist of connected elements in a neighborhood of K. Let the patch 
have finite size such that 


diam(Ar(A^)) < Ch^ 


(3.1) 


where C is a mesh independent constant. 

In a triangle mesh a patch Af{K) typically is the standard four element patch illustrated in 
Figure |l(a)| consisting of K and the three elements neighboring K. For elements neighboring 
the boundary the patch will include a ghost element outside the domain for each element edge 
belonging to the boundary, see Figure [T(b)j As defined in Section 2.1.1 the locations of the 
ghost nodes are set such that the ghost elements are anti-symmetric with respect to K, thus 
preserving properties of structured meshes. 


3.2 Morley Reconstruction 

It is well known that the nonconforming Morley element IfTOll shows optimal convergence in 
the approximation of the Kirchhoff-Love plate bending equation. As noted in [[8]| this element 
is naturally derived in the setting of dG methods for the plate equation by letting /3 —oo in 
( |2.30[ ). An advantage of reconstruction using Morley basis functions is that the jump in the 
normal derivative at the edge midpoint per definition is zero which results in that all interior 


and exterior edge terms E G S\Sc disappear in the bilinear form (2.30). 
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r 


(a) Standard four element (b) Patch on boundary, 
patch. 


Figure 1: Standard patches for (a) interior elements and (b) elements neighboring the boundary. 

The Morley basis functions are constructed so that the deflection field is continuous at the 
nodes and the gradient in the normal direction is continuous at each edge midpoint xe- Clearly 
this is a subspace of WV and as such we define the space of Morley functions 


VVM = {v e WV : [v,„] 1^^ = 0 for all £ G ^ \ {£s U } 


(3.2) 


We define the reconstruction of the normal gradient at an element edge to be the average 
normal gradient of the two neighboring linear triangles. Let V(^) be the set of nodes for ele¬ 
ment K and let Ve{K) be the set of edge midpoints for element K. The reconstruction operator 
TZ : CVi ^ VVM is defined by (TZu )\k = (Vku)\k where TZk '■ CV\{J\r{K)) —)■ V 2 {N'{K)) is 
defined as follows 



.r G V(i^) 
xeVE{K) 


(3.3) 


Next, we will show that this choice of reconstruction yields a method equivalent to the Basic 
Plate Triangle. 

3.2.1 Equivalence with Basic Plate Triangle 

The Basic Plate Triangle (BPT) presented in ifT^fTSll is a triangular plate element using contin¬ 
uous piecewise linear deflections and is derived by combining finite element and finite volume 
techniques. We will now describe our interpretation for derivation of the BPT in the presented 
setting, whereafter we will show equivalence with the method produced by the above choice of 
Morley reconstruction. 

The BPT is a mixed interpolation method where the curvatures and moments are approxi¬ 
mated using piecewise constant functions and the deflection field is approximated using func¬ 
tions in CVi. The fundamental idea in this derivation is that by using partial integration of the 
curvatures such that 


{Kij{u)V)K = {u^ijV)K = 1,2 


(3.4) 


and equivalently for the moments 


{Oij{u),l)K= ^,7 = 1,2 


(3.5) 
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these terms ean be estimated using a CV\ deflection field. 


Starting with the element contribution to the bilinear form (2.30) on each element we have 


aK{u,v) = {Oij{u),Kij{v))K 


(3.6) 


By using that the curvatures and moments are assumed constant on the element and applying 
( |3.4| ) and ( |3.5[ ) we get 


aK{u,v) = -^{Oij{u), l)if(tc,-^(v), 1)/^ = ^ {Xu^nSij + ^u^irij, 1)^^ 




(3.7) 


A deflection field U G CVi is then assumed. As the gradient of a continuous piecewise linear 
function is undefined on element edges they are defined as the average gradient of neighboring 
elements 


U,i\E = {U,i), 1=1,2 


(3.8) 


and l ik ewise for the gradient of the test function v G CVi. Note that this definition of the 
gradient on edges makes all edge terms from the bilinear form ( |2.30| ) to be zero in the method, 
except for edges E G Sc, i.e. the clamped boundary. In the derivation of the BPT the normal 
gradients naturally appear due to the partial integration and are thus enforced weakly on the 
clamped boundary. Thus, there is no need to extend patches on clamped edges with ghost 
elements but if we would the boundary condition would read 


(!/,„)= 0, on Ee Sc (3.9) 

The BPT method is formulated as follows: Find U G CV\ such that 

^ dA:(t/,v) =/(v), forallvGCPi (3.10) 

Keic 


where 

MU,v) = ^ (Un) {Ui)nj, 1)^^ (3.11) 

The average gradient of U is constant on each edge which means the integrals are exactly 
evaluated by midpoint quadrature. We get 

iK(U,v) = ,4 ( E + ( E (3-12) 

' ' \EedK J \EedK / 


where xe is the midpoint of each edge. 

We will now show that the proposed method ( |2.35 1 when using the above Morley recon¬ 
struction is equivalent to the BPT ( 3.10[ 3.12). As previously noted reconstructions into Morley 
space give no edge terms in the bilinear form ( |2.30[ ), except for the clamped boundary, so the 
finite element method reads: Find U G CVi such that 


aKi'R'KU ,TZkv) ciEi'R'KU ,TZkv) = 1{'R-kv), forallvGCPi (3.13) 

KeIC Ee£c 
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where is defined in ( |3.6[ ) and asi-r) can be identified in the bilinear form ( |2.30[ ). This 

boundary term allows us to enforce clamped boundary conditions weakly. As (77.t/)^„ = 
in the Morley reconstruction the enforcement of the clamped boundary condition is equivalent 
to p.9[ ) for large enough /3. 

Apart from the difference in how clamped boundary conditions are enforced, there is also 
a difference in how the load is calculated in the two methods. Disregarding this difference 
for now, if we can show that axiU^v) = ax{TlxU,TZxv) for our choice of TZx, the Morley 
reconstruction yields a method equivalent to the BPT. As the reconstructed functions in the 
above equation are quadratic, both curvatures Kij and moments are constant. Thus, we may 
apply the calculations of p.7| ) and yield 


cikO^kU ,TZxv) = — [X{TZxU)^nSij + IJ.{TZxU)^inj,l)gj^[{TZxv)^i,nj) 


JJdK 


(3.14) 


As the gradient (7lxU)j, z = 1,2 is a linear function, the integrals in the expression above are 
also exactly evaluated through midpoint quadrature. Thus, we have 


,TZxv) = I ^2 ^E{^iT^KU)^nSij + lJ.{TZxU)jnj)\^^ j x 

' ' \EedK / 

^hE{{nxv),inj)\\ (3.15) 

EedK ) 

where xe is the midpoint of each edge. Comparing p.l2[ ) with p.l5| ) we see that the methods 
are equivalent if {TZxw)^i\xE — {'^,i)\xE^ z = 1,2 for w e CV i. Looking at the normal component 
of the gradient we have 


{T^Kw)^n\xE = {W,n) \ 


XE 


(3.16) 


by definition of the reconstruction operator TZx- As the reconstructed function TZxw is quadratic 
and equal to w at the triangle nodes we know that the derivative of TZxw at a midpoint xe in the 
tangential direction is equal to the derivative in the tangential direction of the plane defined by 
the triangle nodes. Lfsing that w^t is continuous over element edges we have 

{TlKw),t\xE = w,t\x^ = (wf) (3.17) 

Thus {TIxw)j\xe — {wj)\xe: i which means that the Morley reconstruction yields a 

method equivalent with BPT, apart from the mentioned differences in enforcement of clamped 
boundary conditions and in load calculation. 


3.3 Fully Quadratic Reconstruction 

For this reconstruction operator we consider for each triangle K the neighborhood N'iK) of 
triangles that share an edge with K. Let V(AA(.fir)) be the set of nodes in N'iK). Then we define 
{7lu)\x = {'R,ku)\x where TZx '■ CVi{N‘{K)) —)■ V 2 {.N‘{K)) is defined as follows 

TZx'■ {TZkv){x) = v{x) , x&V{N'{K)) (3.18) 

In general, except for some special configurations of the nodes in N'{K), this is a well posed 
problem. 
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Xa Xd 



Figure 2: Illustration of two neighboring elements on a structured mesh. 


3.3.1 Relation to Morley Reconstruction 

Consider the notation in Figure]^ We define a structured mesh to be a mesh where the midpoint 
between and will be a criterion which we may formulate as 


XE = 


Xa+Xc Xb+Xd 


(3.19) 


A quadratic function with known values at Xa and Xc will at the midpoint xe have a tangential 
gradient equal to the slope of a linear function with the same known values at Xa and Xc- As 
this is valid for the quadratic polynomials associated with both and K the jump in the 
tangential gradient for these polynomials is zero at the midpoint. The same reasoning is true 
for the points Xb and Xd, and since the midpoint is the same on structured meshes, we conclude 
that the jump in the gradient is zero at xe- Thus, for a structured mesh all interior edge terms 
disappear in ( |2.30| ) as midpoint quadrature exactly evaluates these terms. In this case the fully 
quadratic reconstruction is identical to the Morley reconstruction as the gradient at xe is contin¬ 
uous in both cases. Given a structured mesh, any theoretical results based on the fully quadratic 
reconstruction is thus applicable to the Morley reconstruction/BPT-element. 


3.3.2 Degenerate Patch Configurations 

While it is unlikely that quality mesh generation will produce patch configurations where the 
fully quadratic reconstruction fails, we have identified two possible configurations of the stan¬ 
dard patch where the fully quadratic reconstruction does fail. We call these degenerate patch 
configurations. 

It is possible that two elements neighboring K share two nodes as illustrated in Figure [3(aj| 
and thus only have five degrees of freedom. Obviously this is insufficient for reconstructing a 
complete quadratic polynomial. 

The other degenerate patch configuration occurs when the set of nodes in the patch includes 


four nodes positioned on the same straight line as illustrated in Figure 3(b) Along any straight 
line the quadratic polynomial reduces to a one dimensional quadratic polynomial which is fully 
described using only three nodal values. 

In the next section we will suggest a reconstruction operator that allow extending the patch 
in the case of a degenerate configuration of the nodes. 
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Figure 3: Degenerate configurations of the standard patch, (a) Standard patch only containing 
five nodes, (b) Standard patch where four nodes are positioned along a straight line. 

3.4 Least Squares Fully Quadratic Reconstruction 

To deal with the degenerated cases we consider a larger patch of elements in a neighborhood 
of K and define the reconstruction by exact fitting at the nodes of K and least squares fitting at 
the remaining nodes in the patch. Let V(5) be the set of nodes in a set of elements S. Again 
we define {71u)\k = iJlK^)\K where 71k '■ C-7^\{-^{K)) —^ 7^2{-^{K)) is defined as follows 

f {71kv){x) = v{x) , X e V(i^) 

( min V (CRKv)(n) - v(x)f . Vn = V(V(A’))\V(A’) <3-20) 

y Kv, 

The patch of elements N'iK) is in general the four element standard patch and the above 
reconstruction is then identical to the fully quadratic reconstruction. However, if a degenerate 
patch is detected we extend Af{K) one element at a time using elements neighboring Af{K) 
until the patch is no longer degenerate. 


4 A Priori Error Estimates 

We equip + VW with the following energy norm 

IlkllP = + h\\{Mnniv))\\lK\{gp[jSs) ^ (^">)\\dK\SF 

KeK. 

We note that ||| • ||| is indeed a norm on + VVV since if = 0 then v 

must be a piecewise linear function which due to nodal continuity also is continuous. If also 
'n,KelC = 0 then v is globally linear. Finally, for a well posed problem we 

either need Fc 7^ 0 or that there exists no single straight line Fiine such that F 5 C Fiine- In either 
case we get v = 0 . 

Before turning to our main a priori error estimate we formulate a few lemmas that will be 
needed in the proof. 
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Lemma 4.1. The following inequality holds 


|v||p< 


cEi 

KelC 


kllll, 


forallveH^^ + VVV 


where 


\\ is defined by 


I 1112 _ / —2| |2 I I |2 I /2| |2 I /4| |2 

\v\\\K = h |v|i r + |v|2X + /z IvUic + h \v\ak 


Proof. First recall the well known trace inequality 


(4.2) 


(4.3) 


(4.4) 


which is proven by affinely mapping ^ to a reference element K, using the trace inequalty 
||v||^^ < C'||v|||,||v|| j ^ (see flU), and finally mapping back to K. 


Using the triangle inequality on the interior face contributions of ( |4.1[ ) and then using the 
trace inequality Lemma |43] is readily established. □ □ 


In conformance with (|4.3|) we also define the energy norm for a set of elements S 


Ivllll 


El 

KeS 


\K 


(4.5) 


Furthermore, we will also need to approximate functions using quadratic polynomials on 
each patch. Before we introduce and prove the appropriate estimate for this interpolation error, 
recall the Bramble-Hilbert lemma given in [[5ll . 

Lemma 4.2. (Bramble-Hilbert) Let B be a ball in (O such that (O is star-shaped with respect to 
B and such that its radius p > (1/2) Let Q'^^u be the Taylor polynomial of degree m of u 
averaged over B where u G H'^{co). Then 




m,CO 


k = 0,1, 


(4.6) 


where d = diam{G)) and Ja is the chunkiness parameter of(0. 

Remark. The star-shape criterion on (O means that there should exist a ball B E CO such that 
from any point inside B there is a free line of sight to all points on the boundary of CO. Let pmax 
be the supremum of the radius of all such balls in CO. The chunkiness parameter is then defined 
by 

diam((o) 


rco = 


(4.7) 


We are going to apply the Bramble-Hilbert lemma on each patch, i.e. CO = ^{K) fl /C. 
Further we will need that the chunkiness parameter for all patches is limited and therefore we 
introduce the following restriction on the patches: All patches ^{K) fl /C fulfill the star-shape 
criterion and there exists a global constant 7 such that 

7AA(x)n/c < 7 for all G /C (4.8) 

Note that the shape regularity of the mesh is not sufficient to guarantee ( |4.8| ) for standard four 


element patches. However, in most cases where the standard patch does not comply to ( |4.8| ) we 
may add elements to the patch so that it does. Thus, this restriction will typically not introduce 
any constraints on the mesh. 
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Now we turn to the interpolation error when using quadratie polynomials in the energy 


norm ( |4.5[ ) of a pateh J^iK) fl JC and present the following lemma. 

Lemma 4.3. There is a projection operator P 2 ^k '■ 7/^(A/'(^) fl/C) —)■ V 2 {j^{K) fl/C) such that 


W-Pi,Ku\\\j^{^K)r^]C <Ch (^|M|3,Ar(A:)n/c+^l“l4,A/'(A:)n/c 


(4.9) 


for all sufficiently smooth u. 


Proof By the definition of the energy norm of a pateh ( |4.5[ ) and the quasi-uniformity of the 
mesh it suffiees to show that there exists a patch independent constant C such that 




3-k 


^\ 3 ,Af{K)nK: forfc—1,2,3 


(4.10) 


to prove the lemma. As P 2 ,Kki gives zero contribution to fourth order derivatives the term 


Ch^ A/'(A:)nAC ( |4-9| ) may be directly derived from the definition of the energy norm of a 


patch ( |4.5[ ). 

Next we verify that the requirements of Lemma |4.2| are fulfilled. By restrictions on the 
patches there exists a ball B in every patch such that Af{K) fl /C is star-shaped with respect to 
B. We let the projection operator be defined by the Taylor polynomial of degree 3 of w 
averaged over B, i.e. P 2 ,ku = as defined in [|5l. This will be a quadratic polynomial. 

The constant Cm,a in Lemma [47^ only depends on the domain through the chunkiness pa¬ 
rameter 7a,. As CO = ^{K) n /C we have from the restriction on the patches ( |4.8[ ) that there 
exists a global constant 7 such that 7 (b < 7 for all patches. Using this in the proof of Lemma [ 4 ^ 
in [151 we have that 


Cm,Yo) — Cm,"/ (4-11) 

where Cm,-/ is a patch independent constant and we refer the reader to [[5l for details. 

We complete the proof by applying Lemma |4^ together with ( |4.1 1[ ) which gives 

\^~PTKAkj\r{K)nK:^ |m- 2 “L,A/'(x)nx: (4-12) 

<C2,Y(f’ ^\u\-i,M{K)n]C (4.13) 

(4.14) 


where C is a constant independent of the patch and we used ( |3.1| ) and ( |2.15| ) in the last inequal¬ 
ity. □ 


By Sobolev’s inequality pointwise values are well defined for functions in H^{IC) so the 
Lagrange interpolation operator may be used. We extend the standard Lagrange interpolation 
operator to also define values on ghost elements outside the domain such that n : CPftC) —)■ 
CV\{IC U /Cg). As the functions we need to interpolate lack support outside the domain the 
interpolation values at ghost nodes must be defined. For a ghost node xq associated with 
element K G ICv/e define the interpolation value by 


{nv){xG) = {P2 ,kv){xg) + Akv 


(4.15) 
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where A^rv is given by 

A^V = (v - P 2 .KV) |x=;ci + (v - P 2 XV)\x=X 2 - (v - P 2 ,kv) |x=X 3 


(4.16) 


and the numbering of nodes in K is such that X 3 is the mirror-symmetric node to xq- Note that 
Afcv = 0 for V G V 2 {K). 

We shall also need the following inverse estimate proved in [[8l . 


Lemma 4.4. For all v G VW the following estimate hold 


E'llKW 

KeK. 


IdKMEsUSp) ^ 


<Ci ^(a,y(v),%(v))/c 


(4.17) 


KeK. 


where C denote a constant independent of the meshsize h and the parameter /3. 

Finally, we recall the following lemma from [[8] which we will also give proof to. 

Lemma 4.5. Here we collect three basic results on consistency, continuity, and coercivity: 

1. With u the exact solution of the plate equation and TZU the reconstructed dG solution defined 


by (2.35) we have 


a{u — TZU,TZv) = 0 for all v G C'Pi. (4.18) 

2. There is a constant C, which is independent ofh but in general depends on fi, such that 


a(v,w) < C|||v| 


\w\ 


v,weH^ + VVV 


(4.19) 


(4.20) 


3. For P sufficiently large the coercivity estimate 

c|||v||p < a(v,v) veVF’V, 
holds, with a positive constant c independent ofh and fi. 

Proof. 1. This fact is a direct consequence of the fact that the exact solution u satisfies the 


variational statement (2.29). 


2. Using the Cauchy Schwarz inequality on the definition of the bilinear form (2.30) the in¬ 
equality 


a(v,w) < |||v| 


\w\ 


(4.21) 


immediately follows where |||v|||* is defined by 


KeIC 




dK\gF 


(4.22) 


Estimate ( |4.19[ ) follows by showing that the sum is limited by |||v||p which we prove next. 
We begin by noting that the following equalities hold 


Iklll = \\w - PqwWI + \\Pow\\] 

Ww-PqwWI = ^lk,f||| 


(4.23) 

(4.24) 
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for w G "Pi. As [v,„] is a linear function we may apply these equalities to the first term of the 
sum which together with quasi-uniformity yields 


(4.25) 


^ ^\\b’,n\\\dK\{gFU£s) ^ll^o[^,«]ll5A'\(£:fU£:5) 


where the seconds term already exists in the norm. We can decompose v into v-|-v where v eH^ 
and V G WV. Due to the jump terms in ( |4.25| ) the continuous parts of v give no contribution and 
we may thus replace v with v. To the first term we then apply the triangle inequality to remove 
the jump term and can thereby handle each triangle sharing edge E separately. Applying the 
trace inequality ( |4.4[ ) we get 

^llvilll <c(||v„,||2^+fi2|-,„,|?^^+) <C||v,„j|+ <C(cT,;(v),tc,Xv))x+ (4-26) 

where the last inequality comes from that the Lame parameter ji> 0. 

For the second term in the sum of ( |4.22| ) we begin by subtracting the linear interpolant 
;r[v] = 0. Using the triangle inequality, the trace inequality and interpolation theory we have 




{h ^\\v-1lv\\\++h ^\v-7tv\^f.+^ 


;rv+ill <c{h '^||v-:;rv 


E ^ 

<C|v 
< C{Oij{v), Kjj[v))x+ 


\2 

\2,K+ 


(4.27) 

(4.28) 

(4.29) 


and ( |4.19[ ) is established. 

3. We have 

«(Lv) = '^{Oij{v),Kij{v))K 


KeK. 


- Y. 2((W„(v)),[v,,.))£-j3A-‘||/i,[v,„l|||y£^u£,, (4.30) 

E?,E\{£s'J£f) 


Note that 


((A7nn(v)), [V„])£ = ((M„„(v)),/’o[v,„])i 


(4.31) 


since {Mnn{v)) is a constant and [v,„] is a linear function on E. Using this observation, the 
Cauchy Schwarz inequality followed by the standard inequality lab < £a^ + £^^b^, for any 
positive e, and finally the inverse inequality ( |4.17| ) we obtain 

~ 2((M„„(v)), [v^„])£: > 

Ee£\{£sU£F) 

-eC(a/Xv),?f,y(v))/c-e~^/*~^||PoMll5j^\(£:^u£:f) (4.32) 

KeIC 

Given c, with 0 < c < 1, we choose eC = (1 — c)/3 and take j8 > c -t- we obtain the 
coercivity estimate (4.20). □ 

We are now ready to formulate our main a priori error estimate. 


17 










The final publication is available at: link.springer.com 
Numerische Mathematik (2012) 121:6597 
DOI 10.1007/s00211-011-0429-5 


Theorem 4.6. Assume that the reconstruction operator IZ is linear and satisfies the identity 

nKnv = v, \/veV 2 {M{K)nK:), forallKeK (4.33) 

where TZ is the extended Lagrange interpolation operator. Also assume that u E and that 


the patch restriction (4.8 1 is fulfilled. Then the following a priori error estimate holds 


■77.f/||| < C/r (|m| 3+ /r|M|4) 


(4.34) 


where C is a constant independent of h. 


Before presenting the proof of Theorem |4. 6 1 we remark on how the reeonstruetion operators 
presented in Seetion [prelate to the identity ( |4.33| ) in the theorem. 

Remark. By eonstruetion the fully quadratie reeonstruetion and the least squares fully quadratie 
reeonstruetion satisfy the identity ( |4.33| ). As noted in Seetion [3.3.1 [ this implies that on strue- 
tured meshes the Morley reeonstruetion also satisfies the identity. 

On unstruetured meshes however, the Morley reeonstruetion does not satisfy the identity. 
By the definition of the Morley basis funetions the normal gradient at eaeh edge midpoint must 


be exaetly reeonstrueted if the reeonstrueted quadratie polynomial shall satisfy ( |4.33| ). As we 
in the proposed Morley reeonstruetion use a pair of linear elements to reeonstruet the normal 
gradient on eaeh edge midpoint, we do not have to eonsider the eomplete pateh but rather 
only pairs of elements. To reeonstruet the normal gradient of a quadratie polynomial at the 
edge midpoint in general five degrees of freedom are needed. As we in the proposed Morley 
Reeonstruetion only use four degrees of freedom, an element pair, to reeonstruet the normal 
gradient we generally eannot exaetly reeonstruet quadratie polynomials. 

While this remark does not prove that the Morley reeonstruetion does not eonverge on 
unstruetured meshes it may give some understanding of the numerieal results. 


Proof, of Theorem \4.6\ We first note, using the triangle inequality, that 

|||m — 77.f/||| < |||m — 77.;rM||| + \ \\lZnu — 'RU\\\ 


(4.35) 


where it is the extended Lagrange interpolation operator. Using eoereivity ( |4.20| ), eonsisteney 
( |4.18| ), and the eontinuity properties in Lemma [43] we ean estimate the seeond term as follows 


c\\\TZnu — TW\\\^ < a{'RKu — TZU ,1Zku — TIU) 
= a{TZ7iu — u + u — TZU, IZnu - 
= a{TZ7iu — u, IZnu — IZU) 

< C|||77.;rM — m||| \\\'Rnu — lZU\ 


nu) 


and thus we arrive at 


|||77.;rM — 77.f/||| < C|||m —77.;rM| 


(4.36) 

(4.37) 

(4.38) 

(4.39) 

(4.40) 


Note that the above derivation follows the proof of Cea’s lemma but uses the reeonstruetions 
of the analytieal and finite element solutions, IZltu and TZU, instead of the pure analytieal and 
finite element solutions, u and U. Combining (4.35) and ([4.40) we obtain 


\u — TlU\\\^ < C|||m —77.;rM|||^ < C E |||m — TT-ArTTwIII 

KeJC 


(4.41) 
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where we used Lemma 4.1 in the last inequality. Adding and subtraeting P 2 ^ku and TZk'JiP 2 ^ku 


and then using the triangle inequality we obtain 


\u — TZk'JIu\\\k < \\\u—P2^ku\\\k 

+ \\\P2,KU — TZk'JIP2,Ku\\\k 
-^\\\'T^K'^{P2,KU — u)\\\k 

= /+//+/// 


(4.42) 

(4.43) 


We now eontinue with estimates of Terms I to III. 


Term I. Employing Lemma 4.3 we have 


1= \\\u-P2,ku\\\k < \\\u-P2.ku\\\j^(k)(MC < C(/z|M|3^_^(^)n/c + ^^|M|4,A/'(^)n/c) 


(4.44) 


Term II. Using the assumption (4.33|) on the reeonstruetion operator we eonelude that 


It = \\\P2,KU — TZk'JIP2,Ku\\\k =0 


Term III. Using the following two estimates 

1117r(v - P 2 ,KV) 111 Ar(A-) < C| 11 (v - P 2 ,KV) \ \\M{K)r\K: 

whieh we prove below, we may estimate Term III as follows 


forallvGCPi(Ar(i^)) 

forallve//^(Af(if)n/C) 


III = \\\Rk71{u — P2,Ku)\\\k 

< C\\\n{u - P 2 ,Ku)\\\f^(^K) 

< C|||M-/’2.^:M|||A/-(A')n/c 

< C(/r|M|3^_/v'(/c)n/C + ^^l“l4,A/'(A')n/c) 


(4.45) 


(4.46) 

(4.47) 


(4.48) 

(4.49) 

(4.50) 

(4.51) 


where we used Le mma 4.3 in the last inequality. 

Proof of Estimate (4.46). Let F : ^{K) —>■ ff{K) be a bijeetive eontinuous pieeewise affine 


mapping from a referenee pateh M{K) to the pateh M{K). We note that, due to shape regularity, 
we only need to eonsider a finite number of referenee patehes eorresponding to the different 
topologieal arrangements of the triangles in the pateh. The mapping F takes the form 


Fx = A^x + b^, xEK 




(4.52) 


As F maps a triangle of fixed size from a referenee pateh onto K we have that | detA^| = 
Chi, \\A-^\\<hK/{2p^) < Chfc and by shape regularity ||A^^ || < h^K'^pK) < Chj^^. Next we 

define a mapping IF : CV\{N'{K)) — CVi{N'{K)) by 


V = IFv = voF~ 


(4.53) 
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Together with ( |2.15[ ) we have the estimates 


H„j<C|detA^| '/-||Aj.|r|v|„,(c < CA" 


m,K 
K 


Using (4.54) and (4.551 we conclude that there are constants c and C such that 

. \l<cm 


|V| 

vl,l|~<Clllvl 


where 


iiiV) 


(4.54) 

(4.55) 

(4.56) 

(4.57) 

(4.58) 


Returning to the proof of (4.46) we first show that the inequality holds on the reference 
neighborhood 

~ (4.59) 


III'RkvIIIj^CIIHII^ VveCT>,(A(()C)) 

We note that = 0 if and only if v is constant onN‘{K) but then v = is also constant 

on Af{K) and thus TZkv = v is also constant. Therefore we conclude that |||'^a'v|||^ = 0 if 
-rTTpr = 0 and inequality (|4.59|) thus follows from finite dimensionality. Combining (|4.56|), 

jy{K) 


( |4.57| ) and ( |4.59| ) we get 

WJlKvWk < Clll^vlll J < C|||q||_j^ < C|||v||U(,f, (4.60) 

which concludes the proof of estimate ( 4.46| ). 

Proof of Estimate ( |4.47p . Let w = v — Pi.k^- By contruction of the extended Lagrange inter- 
polant ( |4.15[|4T6 ) and mirror symmetry of ghost elements we have 

(4.61) 

Adding and subtracting w, using the triangle inequality and interpolation error estimates we get 

\\\'ttw\\\\n 



(4.62) 

KciN{K)r\K: 



(4.63) 

KcAf{K)nlC 


<C ^ h '^\w\'i,K+\M2,K 

(4.64) 

KcAf{K)nlC 


< ^ A/'(x)nx: 

(4.65) 


and thus estimate (4.47) follows. 
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We have thereby completed the estimates of Terms I to III in ( |4.43[ ). Using these in ( |4.41| ) 
we thus have 


<C^(/ + // + ///)2 

(4.66) 

KetC 


{h\u\2^yf(K)nlC + h^\^\4jZ{K)nlc) 

(4.67) 

KeJC 


KeIC 

(4.68) 


By shape regularity we have that the number of overlaps in the sum will be finite and thereby 
there exists a constant C such that 


nu\\c <C[h\\l + h^\u\X) <C(/r|M|3 + /r>|4)' 


which completes the proof. 


(4.69) 

□ 


We now turn to an estimate of the norm of the error. This is derived using a duality 
argument (Nitsche’s trick). We assume that for all y/ G -\-WV there is a ^ G such that 

a(v,0) = (v,Va), forallvG//^ + T>PV (4.70) 

and that the following stability estimate holds 

||,(|||4<C||V/|| (4.71) 

On smooth domains and convex bounded polygonal domains where the inner angle at each 
comer is less than 126.3° this assumption is true, see 0. 


Theorem 4.7. If the stability estimate ( 4.71 ) holds, then U satisfies 

\\u — lZU\\ < (|m|3+/z|m|4) 


(4.72) 


for sufficiently regular u. The constant C is independent of h but may in general depend on fi. 

Proof. Setting v = = u — IZU, in the dual problem ( 4.70| ) and using consistency ( |4.18| ) to 

subtract the reconstmction TZn^ of Ttcj) we obtain 


\u — 71U\\^ = a{u — 71U,(I>) 

= a{u — IZU, ^ — TZTt(f) 


<c\\\u-nu\ 


-lZKn\ 


(4.73) 

(4.74) 

(4.75) 


where we used continuity (4.19) in the last step. Next using Theorem |4^ and results ( |4.41| ) in 
its proof we have 


\u-nu\\^ <ch^{\u\^+h\u\^){\<f\^+h\<f\^) 


which together with the stability estimate (4.71) concludes the proof. 


(4.76) 

□ 
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In Theorem |4.6| and Theorem |4.7| we have given a priori error estimates for the reeonstructed 
solution TZU in energy norm and in norm. We now turn to showing an a priori error estimate 
for the eontinuous piecewise linear solution U in norm. 


Theorem 4.8. If the stability estimate {4.71 ) holds, then U satisfies 


W — U\\ < Ch^ (\u\n+ \u\n+h\u\ 


(4.77) 


for sufficiently regular u. The constant C is independent ofh but may in general depend on /3. 
Proof Using triangle inequality we have 


u-U\\ < ||M-7^t/|| + ||7^t/-f/|| 


(4.78) 


where the first term is evaluated by Theorem 4.7 
interpolation estimate 


For the second term we use a standard 


nu-u\\ = \\nu -nlZUW <ch^\nu\ 2 <ch^{\u-nu \2 + \u\ 2 ) ( 4 . 79 ) 


where we in the last inequality use the triangle inequality on the seminorm. 

As the Lame parameter /r > 0 there exists a constant C such that 

\u-nu\l = \u-nu\Y<c\\\u-nu\f (4.80) 

KeK. 

which is limited by Theorem |4.6[ This gives the error estimate 

IIm —f/|| <Ch^{\u\2 + \u\2i + h\u\4) (4.81) 


which concludes the proof. 


□ 


5 Numerical results 

Numerical results will be presented for the following proposed methods: Morley reconstruc¬ 
tion, fully quadratic reconstruction, and least squared fully quadratic reconstruction. Also, for 
comparison we will present results for: the Basic Plate Triangle, the nonconforming Morley 
triangle, a quadratic continuous/discontinuous Galerkin method featuring C® continuity, and a 
quadratic discontinuous Galerkin method continuous at the mesh nodes. 

Note that for the reconstruction methods, the pointwise error is defined as e = u — TZU 
unless otherwise stated. For other methods the pointwise error is as usual defined as e = u — U. 

5.1 Model Problems 

To study the convergence properties of the proposed methods we use two model problems 
where analytical solutions are known. 
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(a) Structured mesh. 


(b) Unstructured mesh. 


Figure 4: Two example triangulations of the unit square with comparable mesh size h. The left 
triangulation (a) is a structured mesh and the right triangulation (b) is an unstructured mesh. 


5.1.1 Problem 1: Simply Supported Plate under Sinusoidal Load 

Consider a simply supported unit square plate, Q. = [0,1]^, with D = \ and v = 0. Find the 
deflection u given the sinusoidal load 

/ = 25ti^ ?,m{nx)sin{2ny) (5.1) 

This problem has the analytical solution u = sin(;rx) sin(2;ry). 

5.1.2 Problem 2: Mixed Boundary Conditions with Uniform Load 

Consider a unit square plate, Q. = [0,1]^, with two opposite sides simply supported, one side 
clamped, and the last side free. Given E = 10^, t = 0.01, v = 0.3 and a uniform load / = 1, 
find the deflection u of the plate. An analytical solution in the form of a series expansion is 
given in Example 46 in ffT^ . 


5.2 Mesh 


The triangulations we consider include both structured and unstructured meshes. The structured 
meshes conform to the criteria discussed in Section 3.3. 1[ Example triangulations of the unit 
square for both structured and unstructured meshes are illustrated in Eigure|^ 


5.3 Numerical Examples 

To illustrate interesting features of the proposed methods we here give a few numerical solu¬ 
tions. 
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Figure 5: Example dG solution to Problem 1 using fully quadratic reconstruction on a coarse 
unstructured mesh. 


5.3.1 Nodal Continuity and Continuity of Normal Gradient 

A reconstructed solution to Problem 1 on a coarse mesh is presented in Figure Note that 
continuity of the nodes is strongly enforced and the continuity of the normal gradients on edge 


midpoints is weakly enforced through the dG method’s ( |2.30| ) inherent penalization of jumps 
in the normal gradient. 


5.3.2 Solution on Mesh including Degenerate Patch 

To illustrate the need of the least squared fully quadratic reconstruction we use a mesh which 
include a degenerate patch, see Figure [6(ajj This mesh was modified to include this patch as it 
is unlikely that degenerate patches appear when using quality mesh generation. The collapsed 
solution is shown in Figure [7^ By extending the patch as in Figure |6(bj] the least squares 
fully quadratic reconstruction gives an accurate solution, see Figure [7(b)| 


5.4 Convergence 


We consider convergence in both the energy norm ( |4.1[ ) and in the I? norm. As the noncon¬ 
forming Morley plate can be viewed as a special case of the quadratic discontinuous Galerkin 


method continuous at the nodes where /3 —)■ oo in ( |2.30| ), the energy norm is applicable also to 
this element. 


5.4.1 Comparison of Morley Reconstruction and Basic Plate Triangle 


As shown in Section 3.2.1 the major difference between the Morley reconstruction and the 
Basic Plate Triangle W2\ lies in the calculation of the load vector. A comparison of the two 
methods using Problem 1 on a structured mesh is shown in Figure and clearly indicate a 
better convergence rate when using the load calculation of the reconstructed Morley method. 
The difference in enforcement of clamped boundary conditions does not produce any noticible 


24 




















The final publication is available at: link.springer.com 
Numerische Mathematik (2012) 121:6597 
DOI 10.1007/s00211-011-0429-5 



(a) Degenerate four triangle patch. (b) Extended patch. 


Figure 6: Example mesh that includes a degenerate patch indicated in (a) and an extension of 
that patch indicated in (b). 


difference in numerical results. While keeping the difference in convergence rate in mind, we 
will from here on let the results for the Morley reconstruction method also represent the beviour 
of the Basic Plate Triangle. 


5.4.2 Convergence on Structured and Unstructured Meshes 


As noted in Section 3.3.1 the Morley reconstruction and the fully quadratic reconstruction 
coincide on structured meshes, and should thereby produce identical results. This is seen in the 
convergence plots for structured meshes. Figures |^and 10, where their paths overlap. 

On unstructured meshes the Morley reconstruction/Basic Plate Triangle does not converge 
to the analytical solution. This is seen in Figures [TT]fT4} As noted in Remark the Morley 
reconstruction does not fulfill the assumption of Theorem |4.6| on unstructured meshes, and thus 
the a priori error estimates are not valid. On the other hand, the fully quadratic reconstruction 
does show optimal convergence on unstructured meshes, as predicted by the a priori estimates. 
In the figures slopes close to 1 for the error in energy norm and slopes close to 2 for the error 
in I? norm indicate optimal convergence. With the noted exception of the Morley reconstruc¬ 
tion/Basic Plate Triangle on unstructured meshes. Figures [9p4| indicate optimal convergence 
for all the compared methods. 

We have previously mentioned that the nonconforming Morley triangle can be seen as a 
special case of the quadratic nodal continuous discontinuous Galerkin method. This is natural 
as the /3 penalty parameter in the dG method enforces continuity of the normal derivatives 
over each edge midpoint, which is the very definition of the Morley basis functions. As shown 
in Figures [9p4l the convergence results for the respective method are close to identical for 
/3 = 100 as used in these calculations. 
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X 10'® 



(a) Collapsed solution due to degenerate patch. 



(b) Accurate solution using extended patch. 

Figure 7: Numerical solution for a simply supported plate under uniform load using LSFQ- 
reconstruction including the degenerate patch is shown in (a) and including the extended patch 
is shown in (b). 
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Figure 8: The error in the numerical solution of Problem 1 versus the mesh size h on structured 
meshes. Slopes for the Basic Plate Triangle and the Reconstructed Morley are 1.28 and 1.99 
respectively. The error e = u — U is measured in norm. 



Figure 9: The error in the numerical solution of Problem 1 versus the mesh size h. Structured 
meshes are used and the error e is measured in energy norm. Note that the reconstructed Morley 
and fully quadratic reconstruction produce identical results. 
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log{h) 


Figure 10: The error in the numerical solution of Problem 1 versus the mesh size h. Structured 
meshes are used and the error e is measured in L? norm. Note that the reconstructed Morley 
and fully quadratic reconstruction produce identical results. 



log(h) 


Figure 11: The error in the numerical solution of Problem 1 versus the mesh size h. Unstruc¬ 
tured meshes are used and the error e is measured in energy norm. 
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log{h) 


Figure 12: The error in the numerical solution of Problem 1 versus the mesh size h. Unstruc¬ 
tured meshes are used and the error e is measured in norm. 



log(h) 


Figure 13: The error in the numerical solution of Problem 2 versus the mesh size h. Unstruc¬ 
tured meshes are used and the error e is measured in energy norm. 
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Figure 14: The error in the numerical solution of Problem 2 versus the mesh size h. Unstruc¬ 
tured meshes are used and the error e is measured in norm. 


5.4.3 Number of Degrees of Freedom 

To give some indication of the performance of these elements in regards to how many degrees 
of freedom are needed to represent the solution we give Figure While it is seen in Figures 
|9 14 that the quadratic cG/dG method has the best performance among the tested methods with 
respect to mesh discretization, Figure[^indicates that the fully quadratic reconstruction has the 
most compact representation performance wise. This is natural as we have smooth solutions. 
Even though the quadratic nodal continuous dG method produce results close to identical to 
those of the Morley triangle with regards to mesh discretization it does feature two degrees of 
freedom on each edge midpoint compared to one for the Morley triangle, explaining that more 
degrees of freedom are needed for par performance. 


5.4.4 Size of penalty parameter j8 


In Figure 16 we present some numerical results for various /3. As might be suspected the 
fully quadratic reconstruction exhibits locking effects when j8 is to large. This is natural as 
neighbouring elements share much of the information through the patch construction. A more 
surprising result is that the quadratic cG/dG method does not seem to exhibit such locking ef¬ 
fects for large j8. This indicates that the finite element space of continuous piecewise quadratic 
polynomials with continuous normal gradients on edge midpoints is large enough to accurately 
approximate the solution. If we on the other hand change the projection operator in the penalty 
term from the projection onto constants Fq to the projection onto linear functions P\ the cG/dG 
method exhibits locking effects for large /3. 


A mesh independent lower bound for /3 can be calculated if a suitable choice of h in ( |2.30| ) 
on each edge is made, see [Q . However, for the numerical results in this paper we have used a 
global mesh size parameter for h. As the meshes used in the numerical results in this paper are 
quasi-uniform this should be sufficient. 
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Figure 15: The error in the numerical solution of Problem 2 versus the number of degrees of 
freedom needed. Unstructured meshes are used and the error e is measured in energy norm. 



\om 


Figure 16: The error in the numerical solution of Problem 1 versus the mesh size using various 
j8. Solid lines indicate j8 = 10^, dashed lines indicate j8 = 10^, and dash-dot lines indicate 
j8 = 10^. Unstructured meshes are used and the error e is measured in energy norm. 
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